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Abstract 

Quark propagators with arbitrary sources and sinks can be obtained more 
efficiently using a pseudofermion method with a mode-shifted action. Mode- 
shifting solves the problem of critical slowing down (for light quarks) induced 
by low eigenmodes of the Dirac operator. The method allows the full physical 
content of every gauge conhguration to be extracted, and should be especially 
helpful for unquenched QCD calculations. The method can be applied for all the 
conventional quark actions: Wilson, Sheikoleslami-Wohlert, Kogut-Susskind, as 
well as Ginsparg-Wilson compliant overlap actions. The statistical properties 
of the method are examined and examples of physical processes under study 
are presented. 
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1 Introduction 


The development of more powerful computing platforms as well as improvements in 
algorithms suggest that unquenched lattice QCD gauge conhgurations on reasonably 
large lattices should become available in the not too distant future. The generation of 
equilibrated and decorrelated dynamical gauge conhgurations for lighter quark masses 
will nevertheless remain extremely costly in computational terms, and elementary 
considerations of load-balancing suggest that we should be willing to expend a cor¬ 
respondingly large computational effort in extracting the maximum physical content 
from each available gauge conhguration. Hadronic correlators built from conven¬ 
tional quark propagators (which give the quark propagation amplitude from a single 
source to all points on the lattice, for example) evidently exploit only a fraction of 
the physical information latent in each gauge conhguration. In this paper we explore 
the statistical properties and computational feasibility of an alternative stochastic 
approach to generating quark propagators, one which supplies simultaneously quark 
propagation amplitudes from any point on a space-time lattice to any other (i.e. 
all-point propagators). The use of pseudofermion helds has been studied previously 
0, particularly for single light quark systems (e.g. heavy-light mesons). With 
these approaches, as the quark mass becomes light the convergence of pseudofermion 
Monte Carlo calculation suhers critical slowing down. In this paper, we observe that 
for physical correlators (in momentum space) the statistical noise problem resides 
mainly in its low momentum behaviour. We will hnd that the convergence can be 
greatly improved by separating oh the lowest eigenmodes of the Dirac operator and 
treating them exactly while doing the Monte Carlo calculations with a mode-shifted 
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action. 


In Section 2 we review the essential properties of pseudofermion fields which allow 
the computation of an all-point propagator. Although the detailed results presented 
in this paper primarily concern Wilson or clover-improved quark actions, we also 
indicate here how the method can be applied to the computation of an all-point 
overlap operator. Of course, the basic properties of pseudofermion helds will already 
be familiar to practitioners of the hybrid MonteCarlo or the Luescher multibo¬ 
son method Q]. In Section 3 we describe in some detail the use of pseudofermion 
quark propagators in constructing a variety of hadronic correlators (corresponding 
to both two- and three-point hadronic Greens functions), such as disconnected parts 
contributing to form factors or to hairpin amplitudes for isoscalar mesons, etc. In 
Section 4 we describe the statistical properties of hadronic correlators computed using 
pseudofermion all point propagators. The availability of all-point propagators allows 
high-statistics evaluations of the full four-momentum structure of such correlators. 
Detailed autocorrelation studies reveal that the method gives very accurate results 
(with rapidly decorrelating amplitudes with a time constant of a few heat-bath sweeps 
of the pseudofermion helds) for all but the lowest lattice momenta. The critical slow¬ 
ing down in these modes is related to the presence of low eigenmodes of the hermitian 
Wilson-Dirac operator. Projection and shifting of a few low-lying modes turns out 
to be computationally straightforward, and dually in Section 5 we show that such 
mode-shifted pseudofermion simulations allow accurate extraction of all momentum 
components of hadronic correlators. 
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2 Pseudofermion Actions and All-Point Quark 
Propagators 

A pseudofermion action suitable for computing all-point quark propagators is con¬ 
structed from the quadratic form Q defining the basic quark action of the theory: 

‘S'quark ax Q a:x.,by'lpby (1) 

ax.,by 

where a, b are spin-color indices and x, y lattice sites. For the time-being, we shall con¬ 
sider Wilson or clover-improved actions only (later, we return to the case of Ginsparg- 
Wilson compliant overlap actions). Then the operator H = 75 Q is hermitian and 
we can form a positive-dehnite bosonic (“pseudofermion”) action as follows 

Sp! = E (2) 

ax,by 

in terms of a complex bosonic pseudofermion held (p. Averages of products of the 
pseudofermion held with respect to the Gaussian weight dehned in Eq(g) yield inverses 
of H'^ in the usual way: 

<< F{p,p*) » = Z-^ j 

Z = j 

« 0ax0by » = (hf“^)ax,6y (3) 

Note that the Gaussian dependence of the pseudofermion action implies 

« 5pf »= iVcoiiVspinG = 12V (4) 

where V is the lattice volume (=number of lattice sites). This exact result is extremely 
useful in checking for errors in the simulation algorithm and in determining equili¬ 
bration of an ensemble. Of course, the quark propagator of interest in lattice QGD 
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simulations is (up to a trivial 75 factor) H not H This can easily be achieved by 
a single additional “dslash” operation applied to the conjugate pseudofermion field: 

^ax06y (-^ )ax,6y ^b)ax.,by (h) 

Note that separate pseudofermion helds are needed for each quark propagator as av¬ 
erages of products of (say) four bosonic pseudofermion fields will produce contractions 
with the wrong sign relative to the corresponding fermionic 4 quark amplitudes. The 
simulation of averages of the kind found in ill) can be readily accomplished by 
heat-bath updates, due to the simple Gaussian dependence of the action on the helds. 
For either Wilson or clover-improved actions, the dependence of the pseudofermion 
weight on the pseudofermion held at a specihc lattice site x takes the form 

where is a complex spin-color vector assembled from the pseudofermion held at 
nearest and next-to-nearest sites, as well as appropriate gauge-link variables connect¬ 
ing to those sites. The 12x12 matrix ^4 is a multiple of the identity (specihcally 
^4 = 1 -|- 16k^) for Wilson actions, and a more complicated hermitian matrix assem¬ 
bled from clover gauge helds in the SW-improved case. In either case, the heat-bath 
update of 0x is readily managed by completing the square in (^. For the clover 
actions, the matrices Axab can be prediagonalized just once at the start of the simu¬ 
lation, and the resulting stored eigenvalues and eigenvectors used to quickly generate 
(f) updates at each site according to the weight (^. 

The pseudofermion method can be readily generalized to study the all-point over¬ 
lap operator arising from an overlap action satisfying the Ginsparg-Wilson condition 
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[^, 1^. Let H be the hermitian Wilson-Dirac operator with suitably chosen overlap 
mass. The all-point overlap operator for arbitrary bare quark mass is trivially com¬ 
putable once all matrix elements of the nonlocal operator e{H) = H ■ are 

obtained. Using an optimal rational approximation |Q , this nonlocal operator can be 
written 

N ^ 

e{H) ~ (ao + ^ , , )H (7) 

where are real and the are real positive. The number of pole terms N needed 
for a given level of uniform accuracy over the spectrum of H is related in a fairly 
straightforward way to the condition number of H, but typical studies of the overlap 
operator have used 10 < < 100. We shall return to this issue in Section 4, where we 

show that mode-shifting can be used to dramatically improve this condition number 
and reduce the number of poles needed. The needed all-point operator can clearly 
be obtained by an average of N pseudofermion helds, m = 1, N: we begin from 
the positive dehnite action 

N 

^pf,overlap = E (8) 

m=l 

and construct the desired combination of pole terms from a corresponding combina¬ 
tion of pseudofermion helds, averaged relative to the weight (H): 

e{H)a^,by =« E«rn</’ixV^”'^6y >> +aoHa^,by (9) 

m 

The computation of all-point overlap propagators will require a simulation within a 
simulation (analogously to the situation for single source overlap propagators, where 
inversions within an outer conjugate gradient inversion are required, 0). 
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3 Hadronic Correlators from All-Point Pseudo¬ 


fermion Quark Propagators 

In general, the computation of multipoint hadronic correlators involving n quark 
propagators can be reduced to convolutions of n pseudofermion helds, rapidly com¬ 
puted by fast Fourier transform (FFT). In fact, using FFT we can easily construct 
a wide range of correlators involving both local and smeared operators. In this sec¬ 
tion we shall illustrate this by considering a number of examples of physical interest. 
We note here that it may be advisable to combine all-point with conventional hxed 
source (or sink) propagators (obtained, say, from a conjugate gradient inversion), as 
accurate results become increasingly difficult as the number of allpoint propagators 
used increases, due to the large condition number of H, as will become clear in subse¬ 
quent sections. Furthermore, the projection methods described in Section 5 become 
increasingly essential in reducing the statistical errors of the pseudofermion average 
as the number of independent all-point propagators increase. 

3.1 Local Two-point Correlators 

The all-point propagators obtained by the pseudofermion technique can be used to 
construct the complete momentum-space dependence of 2-point correlators of scalar 
or pseudoscalar densities, or vector or axial-vector currents, while exploiting the full 
physical content of each gauge conhguration. For example, we may be interested in 
the full 4-momentum transform 

Ap5-P5(g)=Ee*‘''^"“"A(a;,y) 

x,y 
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of the 2-point pseudoscalar correlator A{x,y), given by 

A{x,y) = < 0 |T{'h(a;) 75 d'(a;) ^(?/) 75 d'(y )}|0 > (10) 

= - < tr((Q“S5)xy(Q“S5)yx) > (H) 

^ H)yijXvb{X^H)xa 

ab 

~ {4^^Hx)yy{X^H(j))xx >> 

~ (0 ■ X)yy{X ■ 4’)xx >> (12) 

where (10) represents a conventional operator VEV, the single brackets in (11) refer to 
a functional average over gauge-helds, and the double brackets in (12) imply averages 
over gauge conhgurations as well as a pseudofermion average for each gauge conhgu- 
ration to determine the all-point propagators Q~^ in (11). The full momentum-space 
correlator thus becomes an easily evaluated fast Fourier transform of products of 
pseudofermion helds: 

Aps-ps{q) = - « FFT(x • 0)(q)FFT((^ • x)(-q) >> (13) 

At this point, it may be of use to the reader to indicate the computational require- 
mentsof these pseudofermion simulations (see [Ql for more details). For a lattice, 
a single heat-bath update of the two pseudofermion helds 0, x requires 0.366 sec. 
on a 1.5 GHz Pentium-4 processor. The convolutions and FFT operations required 
to obtain the desired four-momentum held Aps-ps{q) in ( P!3[ ) require an additional 
0.024 sec. and are performed after every 2 heat-bath updates of 0, y. Typically, a 
sufficiently accurate pseudofermion average for Aps-ps{q) was obtained from 20000 
measurements, corresponding to 2.1 Pentium-4 hrs. For 10^x20 lattices, a heat-bath 
update costs 3.30 sec, a measurement of A(g^) requires 0.43 sec. (again, performed 


every 2 heat-bath steps), and hnal averages are taken from 7000 measurements, cor¬ 
responding to 6.8 Pentium-4 hrs. For comparison, the evaluation of a conventional 
conjugate-gradient single-source propagator on a 10^x20 lattice requires 1.2 Pentium-4 
hrs. 


3.2 Smeared Meson Propagators 

In lattice QCD, hadron spectroscopy is usually studied using smeared hadron 
sources/sinks to optimize the ground-state signal in each hadronic channel. A typical 
multistate propagation amplitude from Euclidean time ti to tj might therefore involve 
evaluation of the matrix 


= (14) 

Ui^)fl3{y) <0\T{^{z + X,ti)'y 50 a'^{z,ti)^{w + y,tf)j 50 f 3 ^{w,tf)\ 0 > 

xyzw 

where fa{,x) are a set of spatial smearing wavefunctions, 'k(x, f) denotes the quark 
held at spatial (lattice) point x and time t, and Oa are the appropriate spin matrices 
for the desired hadronic channels. The sums over z, w project the physical states onto 
zero-momentum. Of course, delta function choices for fa allow us to use local sources 
or sinks. In terms of pseudofermion quark propagators 


■Mafsitijtf) — fa{^)ff3{y){^^fw+y,tfOf}Xw,tfXz+x,tiOa(j)z,ti >> 

xyzw 

<< fz+x,tiOafz,ti\w+y,tfOpXw,tf ( 1 ^) 


where the second line of (|^ (77j=number of havors) is present only for isoscalar 
amplitudes where there is a nonvanishing loop-back (“hairpin”) amplitude. The four¬ 
fold spatial sum in (pA]) can fortunately be reduced to operations linear in the spatial 
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volume Vs of the lattice via the magic of Fourier transforms. Dehne smeared pseudo¬ 
fermion helds 




s ^ 

S p 


where /a(p) is the Fourier transform smearing function and <pp^t represents the spa¬ 
tially Fourier transformed (and time-sliced) pseudofermion held. Taking the isovector 
amplitude part of (|^ for simplicity, one hnds that this reduces to 


A A isovec 






» 


( 16 ) 


The isoscalar contribution, if present, is trivially obtained by interchanging the 
psendofermion helds in an obvious way. 


3.3 Three Point Functions - The Pion Form Factor 

The computation of more complicated correlators, snch as the 3-point functions 
needed to extract form factors, is greatly facilitated if we have all-point qnark propa¬ 
gators at our disposal. A typical example is shown in Fig.l, where a smeared meson 
created at time to propagates to time ti, where spacelike momentum q is injected 
at spacetime point {y, ti) by an electromagnetic current, followed by propagation of 
the hnal-state meson to time t 2 , when it is removed by an appropriate smeared-sink 
operator. 

The quark diagram displayed in Fig. 1 represents the connected contribution to 
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Figure 1: Quark flow diagram for computing pion form factor(connected part) 
the following hadronic correlator: 

j.,.,.,(<!) = E 

wxyz 

where in this case the same smearing function (optimized for a pion, say) is 
applied at both initial and hnal times. The complete momentum-dependence of this 
3-point function can be obtained by using an all-point pseudofermion propagator 
for the quark propagator from (y, fi) to {x,t 2 )- The quark propagators into the 
point sink at (0,to) and from the smeared source point at (tc, to) can be computed 
by conventional conjugate-gradient techniqnes. Again, by using Fourier transforms 
appropriately, the calculation of Jtotit 2 {^ can be rednced to operations at most linear 
in the spatial lattice volume. Finally, as for the case of smeared meson propagators 
discussed above, the disconnected contribution involving the contraction of 
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back to h) can be calculated with no extra effort, as the all-point propagator for 
just this case has been computed. In this case the pseudofermion held will appear in 
an average of the form << (f>ffti 7 o 754 ’yti »■ 
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4 Statistical Properties of Pseudofermion Propa¬ 


gators 

In this section we shall describe the results of detailed studies of the statistical prop¬ 
erties of pseudofermion propagators computed from the basic formulas (3). We hrst 
consider a case where only one all-point propagator is used in assembling the full 
hadronic amplitude- namely, the 3-point correlator J{q) giving the pion formfactor 
discussed in Section 3.3. We shall see shortly that as a consequence of the typically 
high condition number of the hermitian Wilson-Dirac operator H, autocorrelations 
for low-momentum amplitudes become progressively longer as more all-point propa¬ 
gators are introduced into the calculation (in Section 5, we shall show how to £x this 
problem by mode-shifting). Here, we begin with an application where autocorrela¬ 
tions are relatively unproblematic. 

For the pion form-factor calculation, we generated quenched configurations at 
f3 =5.9 on 12^x24 lattices, and studied the resulting quark propagators at k =0.1590 
(with the Wilson action). The results described here were obtained by studying the 
simulation of the quantity J(g) for a randomly chosen gauge configuration from this 
ensemble: examination of other configurations reveals that the behavior we describe 
is generic. From Eq(4) the (infinite ensemble) average pseudofermion action should 
therefore be Fpf =497664. From Fig.2 it is apparent that the equilibrium value for 
this quantity is reached rather quickly, after about 20 heat-bath sweeps through the 
lattice. An average of the value of 200 consecutive pseudofermion configurations after 
this gives for example Fpf =497599 with a standard deviation of ~650. Typically, we 
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Figure 2: Equilibration of pseudofermion average on 12^x24 lattices 

have performed measurements using pseudofermion fields after 100 initial sweeps. 

The decorrelation of hadronic amplitudes in the course of a MonteCarlo simu¬ 
lation of the pseudofermion propagator (|^) is extremely sensitive to the particular 
momentum component being calculated. In particular, low momentum components 
have a large overlap with the smallest eigenmodes of H, which can have very small 
eigenvalues (~10“^ is not uncommon). These low modes decorrelate only slowly in 
any local update of the action (2). This property becomes immediately apparent 
when we examine the momentum dependence of either the convergence of cumulative 
averages (Fig.3) of J{q^) (taking t 2 — h = ti — to=3) or the autocorrelation function 
of the same quantity (Fig.4), as a function of number of pseudofermion heat-bath 
sweeps performed. The autocorrelation time (defined as the integral under the auto¬ 
correlation curve of Fig. 4) turns out to be about 13 for the zero-momentum mode 
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Figure 3: Cumulative averages for pion form-factor amplitude = 0,1 

and 7 for the =1 (lattice units) mode: higher momenta (not shown here) yield 
autocorrelation times of order unity. In fact, the autocorrelation times for this cor¬ 
relator are fairly mild in comparison to the cases involving two all-point propagators 
discussed below: 8000 pseudofermion sweeps (= ~ 12 hrs on a 1.5GHz Pentium 4) 
allow J(g) to be extracted with error bars well below the intrinsic fluctuation of the 
correlator from one gauge conhguration to the next (see Fig.5). 

The critical slowing down seen in low eigenmodes becomes a more serious problem 
in situations where two separate all-point propagators are required to construct the 
desired correlator, as in the local two-point correlators discussed in Section 3.1. In 
this case the very large condition number of implies an extremely slow convergence 
for the zero-momentum component of the correlator. For example, from (11) we see 
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Autocorrelation of J(q^) 







that the local pseudoscalar density correlator can be written, in momentum space 

Aps-ps(«) = E (18) 

which becomes, for zero momentum 

Ap5-ps(0) = Tr(^) (19) 

As one not uncommonly encounters small eigenvalues of H, it is apparent that a 
few low eigenmodes can contribute disproportionately to this quantity. Moreover, 
these are exactly the modes that decorrelate most slowly in the pseudofermion sim¬ 
ulation. To illustrate this, we have studied P] hadronic 2-point correlators on an 
ensemble of unquenched configurations generated with the truncated determinant al¬ 
gorithm (TDA,[^]) on physically large coarse 6"^ lattices (lattice spacing ~0.4 fm). 
Cumulative averages for Aps-psiq) for a range of values of for a typical gauge 
conhguration in this ensemble, are shown in Fig. 5. The zero momentum mode is 
clearly not convergent even after 8000 pseudofermion sweeps, while even the smallest 
nonzero (lattice) momentum component shows much more rapid convergence. For 
the particular gauge conhguration illustrated here, the lowest eigenvalue of H turns 
out to be 0.0024, which contributes 54% of the total zero-momentum value Tr-^ 

! The problem for this lowest mode can be seen (Fig. 6) in another guise in the 
autocorrelation curves for Aps-ps{q^),0 < < 2 (the autocorrelation time is ~ 1 

pseudofermion sweep for > 2 (lattice units) so these curves are not shown). For¬ 
tunately, the critical slowing down experienced in these pseudofermion simulations of 
all-point propagators appears only to infect the very lowest momentum components. 
We shall see in the next section that the problem can be eliminated for these com¬ 
ponents by mode-shifting a relatively small number of low eigenmodes of H. This 
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Figure 5: Cumulative averages of pseudoscalar density correlator Aps- psiq"^), 0 < 
< 4 

preconditioning substantially reduces autocorrelation times and allows us to extract 
reasonably accurate values even for the zero-momentum component of hadronic cor¬ 
relators. The extraction and separate treatment of low eigenmodes is also essential 
in calculating accurate all-point overlap operators. 
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5 Mode-Shifted Simulations of All-Point Propaga¬ 


tors 


The slow convergence of the lowest momentum modes discussed in the preceding sec¬ 
tion can be substantially accelerated by shifting the low eigenmodes of the Hermitian 
Wilson-Dirac operator H responsible for the critical slowing down. In the case of 
the two-point correlator (pISD , the relevant parameter is the condition number of H^, 
which (for the ensemble of unquenched 6^ lattices discussed in Section 4) can typically 
be reduced by two orders of magnitude by shifting the lowest 10 eigenmodes of H^. 
More generally, define 

N 

Hs = H + Y,6iwM ( 20 ) 

i=l 

where Vj are a complete orthonormal set of eigenmodes of H, Hvi = AjVj, and the 
lowest N modes (in absolute value) are shifted: 

= aS*) - A, (21) 


For simplicity we shall take = sign(Aj-^^) henceforth, although any value with 
magnitude of order unity will do. The extraction of low eigenmodes of H is computa¬ 
tionally straightforward using implicitly restarted Arnold! techniques each mode 
requires a few minutes on a Pentium-4 processor for the 6^ lattices discussed here. 


Corresponding to the shifted operator Hg defined in (|20|) is a shifted pseudofermion 
action: 


S.,pt = E €MU.i.y<f».y (22) 

ax,by 

Once the relevant low eigenmodes Vi are known, the heat-bath update of the pseudo¬ 
fermion field (f) can be performed by trivial modifications of the procedure outlined 
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in Section 2. The added computational load is not large: if iV=10 modes are shifted, 
the time required for a pseudofermion update increases by about 20%. For the rest 
of this section, we shall use the double-bracket notation << .... >> introduced in 
Section 2 to indicate averages relative to the weight generated by the shifted action 
Ss,pi- Accordingly, the unshifted quark propagator is given by 

N 

^ax,by faxfby >> ^i'^i,ax^i,by (23) 

i=l 

where 

A. = 1 - ^ (24) 

and the tilde notation on the right-hand-side of (p3D still refers to the unshifted 
operator H, as in (|^). The pseudoscalar correlator in (10-12) can therefore be written 

Aps_Ps(9) = E 

ax,6y 

= { « (0 ■ x)y{x ■f)x-Y ■ Vi)x(Vi • x)y + (Vi • (j))x{f ■ Vi)y) » 

x,y i 

+ ^ AiA^-(vj ■ Vi)x(v, ■ Vj)y} (25) 

id 

The term involving a double sum over shifted modes in (|^) does not involve 
pseudofermion helds and is therefore calculated just once. We see that the usual result 
(12) has to be supplemented by an average of overlaps of the two pseudofermion helds 
with the shifted eigenmodes. Again, this is computationally perfectly manageable. 

The dramatic effect of shifting even a few low eigenmodes on the convergence of 
correlators built from pseudofermion averages is illustrated in Figs. 7,8,9,10. As ex¬ 
pected, the worst behavior is found in the zero momentum mode, where the unshifted 
averages (Fig. 7) are still 60% below the correct answer after 50,000 pseudofermion 
sweeps, while only 10,000 sweeps already give a reasonably good result after the 10 
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Figure 7: Cumulative average of zero-momentum pseudoscalar correlator 

lowest modes are shifted. The autocorrelation curves for the zero-momentum compo¬ 
nent tell the same story (Fig. 8): the autocorrelation time is about 140 sweeps for the 
unshifted simulation, dropping to 11 (5) sweeps after shifting 10 (resp. 20) modes. 
The overall situation is much better for the lowest non-zero momentum mode, = 1, 
as shown in Figs. 9,10. Here the autocorrelation times are roughly 8,3,2 sweeps for 
simulations with 0,10 and 20 shifted modes respectively, and the cumulative averages 
reflect a correspondingly higher stability. The numerical evidence from these simu¬ 
lations clearly suggests that mode-shifting with 77=10 modes, for the ensemble of 6^ 
TDA lattices considered here, is perfectly adequate for obtaining accurate results at 
all momenta. 

The calculation of an all-point overlap operator using ( 0 ) can be greatly facilitated 
by mode-shifting. The number N of poles needed in the optimal rational approxi- 
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gure 8: Autocorrelation curve for zero-momentum pseudoscalar correlator 



Figure 9: Cumulative average for unit-momentum pseudoscalar correlator 
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Figure 10: Autocorrelation curve for unit-momentum pseudoscalar correlator 

mation to achieve a desired uniform accuracy for e{H) over the full spectrum of H 
is directly related to the condition number (ratio of highest to lowest eigenvalues) of 
H. On the other hand, provided the mode shifting preserves the algebraic sign of 
the shifted eigenvalues (sign(A^^^)=sign(Ai)), one clearly has e{Hs) = e{H). Accord¬ 
ingly, the condition number of H can be drastically reduced by computing the overlap 
pseudofermion operator using (||) with H replaced with Hg as the pseudofermion ac¬ 
tion: the required code is identical to that used in the mode-shifted calculations of 
conventional Wilson propagators described above. 
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